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Abstract 

The spatially homogeneous, isotropic Standard Cosmological Model 
appears to describe our Universe reasonably well. However, Einstein's 
equations allow a much larger class of cosmological solutions. Theo- 
rems originally due to Penrose and Hawking predict that all such mod- 
els (assuming reasonable matter properties) will have an initial singu- 
larity. The nature of this singularity in generic cosmologies remains a 
major open question in general relativity. Spatially homogeneous but 
possibly anisotropic cosmologies have two types of singularities: (1) ve- 
locity dominated — (reversing the time direction) the universe evolves to 
the singularity with fixed anisotropic collapse rates ; (2) Mixmaster — the 
anisotropic collapse rates change in a deterministically chaotic way. Much 
less is known about spatially inhomogeneous universes. Belinskii, Khalat- 
nikov, and Lifshitz (BKL) claimed long ago that a generic universe would 
evolve toward the singularity as a different Mixmaster universe at each 
spatial point. We shall report on the results of a program to test the 
BKL conjecture numerically. Results include a new algorithm to evolve 
homogeneous Mixmaster models, demonstration of velocity dominance 
and understanding of evolution toward velocity dominance in the plane 
symmetric Gowdy universes (spatial dependence in one direction), demon- 
stration of velocity dominance in polarized U(l) symmetric cosmologies 
(spatial dependence in two directions) , and exploration of departures from 
velocity dominance in generic U(l) universes. 



1 Introduction 

We shall describe a series of numerical studies of the nature of singularities in 
cosmological models. Since the interiors of black holes can be described locally 
as cosmological models, it is possible that our methods and results may be useful 
to the participants in this conference. 
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The generic singularity in spatially homogeneous cosmologies is reasonably 
well understood. The approach to it asymptotically falls into two classes. The 
first, called asymptotically velocity term dominated (AVTD) pi 0], refers to 
a cosmology that approaches the Kasner (vacuum, Bianchi I) solution as 
t — ► oo . (Spatially homogeneous universes can be described as a sequence of 
homogeneous spaces labeled by r. Here we shall choose r so that t = oo co- 
incides with the singularity.) An example of such a solution is the vacuum 
Bianchi II model Q which begins with a fixed set of Kasner-like anisotropic 
expansion rates, and, possibly, makes one change of the rates in a prescribed 
way (Mixmaster-like bounce) and then continues to t = oo as a fixed Kasner 
solution. In contrast are the homogeneous cosmologies which display Mixmas- 
ter dynamics such as vacuum Bianchi VIII and IX and Bianchi VI 

and Bianchi I with a magnetic field J§|, ||, [l(| . Jantzen |0| has discussed other 
examples. Mixmaster dynamics describes an approach to the singularity which 
is a sequence of Kasner epochs with a prescription, originally due to Belinskii, 
Khalatnikov, and Lifshitz (BKL) ||, for relating one Kasner epoch to the next. 
Some of the Mixmaster bounces (era changes) display sensitivity to initial con- 
ditions one usually associates with chaos and in fact Mixmaster dynamics is 



chaotic 12 1 . The vacuum Bianchi I (Kasner) solution is distinguished from the 
other Bianchi types in that the spatial scalar curvature 3 R, (proporional to) the 
minisuperspace (MSS) potential || [l|, vanishes identically. But 3 R arises in 
other Bianchi types due to spatial dependence of the metric in a coordinate ba- 
sis. Thus an AVTD singularity is also characterized as a regime in which terms 
containing or arising from spatial derivatives no longer influence the dynam- 
ics. This means that the Mixmaster models do not have an AVTD singularity 
since the influence of the spatial derivatives (through the MSS potential) never 
disappears — there is no last bounce. 

In the late 1960's, BKL claimed to show that singularities in generic solutions 
to Einstein's equations are locally of the Mixmaster type H. This means that 
each point of a spatially inhomogeneous universe could collapse to the singularity 
as the Mixmaster sequence of Kasner models. (It has been argued that this could 
generate a fractal spatial structure JlB], |l(| |l7).) In contrast, each point of a 
cosmology with an AVTD singularity evolves asymptotically as a fixed Kasner 
model. Although the BKL result is controversial pH , it provides a hypothesis 
for testing. Our ultimate objective is to test the BKL conjecture numerically. 



2 Numerical Methods 

The work reported here was performed by using symplectic ODE and PDE 
solvers While other numerical methods may be used to solve Einstein's 

equations for the models discussed here, symplectic methods have proved ex- 
tremely advantageous for Mixmaster models and have also worked quite well in 
the Gowdy plane symmetric and polarized U(l) symmetric cosmologies po| , pl| , 
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^2, 23 . Consider a system with one degree of freedom described by q(t) and its 
canonically conjugate momentum p(t) with a Hamiltonian 

H=^-+V(q) = H K + H v . (1) 
2m 

Note that the subhamiltonians Hk and Hy separately yield equations of motion 
which are exactly solvable no matter the form of V. Variation of Hk yields 
q = p/m, p = with solution 

p(t + At)=p(t) , q(t + At) = q(t) + — At. (2) 

m 



Variation of Hy yields q — 0, p — —dV/dq with solution 

v 

At. (3) 



q(t + At) = q(t) , pit + At) = p(t) - ^ 

aq 



Note that the absence of momenta in Hy makes @ exact for any V(q). One can 
then demonstrate that to evolve from t to t + At an evolution operator U{%\ (At) 
can be constructed from the evolution sub-operators UjiiAt) and Uy(At) ob- 
tained from @ and (||) . One can show that |lS[ ] 

U (2) {At)=U K {At/2)Uv{At)U K (At/2) (4) 

reproduces the true evolution operator through order (At) 2 . Suzuki has devel- 
oped a prescription to represent the full evolution operator to arbitrary order 
[ p4| . For example 

U w ( At) = U m (sAt) U m [(1 - 2s) At] U m (a At) (5) 

where s = 1/(2 — 2 1 / 3 ). The advantage of Suzuki's approach is that one only 
needs to construct U{2) explicitly. I4(2 n ) is then constructed from appropriate 
combinations of U( 2n -2)- 

The generalization of this method to N degrees of freedom and to fields 
is straightforward. In the latter case, V^[<f(t)] — > V[q(x, t)] so that dV/dq be- 
comes the functional derivative 6V/5q. On the computational spatial lattice, 
the derivatives that are obtained in the expression for the functional derivative 
must be represented in differenced form. We note that, to preserve nth or- 
der accuracy in time, nth order accurate spatial differencing is required. Some 
discussion of this has been given elsewhere . 

3 Application to Mixmaster Dynamics 

(Diagonal) Bianchi Class A cosmologies |l3| are described by the metric 

ds 2 = -e m dt 2 + {e^). 3 ,<tV (6) 
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where /3y = diag(- 2/3+, (3+ + y/3/3-, (3+ - a/3/3-), cfcr 1 = Cj fc er J ' A a k defines 
the Bianchi type and t is the BKL time coordinate. Einstein's equations can be 
obtained by variation of the Hamiltonian 

2H = -p 2 n +pl+ P 2 _ + V({3+,p-,n). (7) 

The logarithmic anisotropic scale factors a, £, and 7 are given in terms of the 
logarithmic volume f2 and anisotropic shears (3± as 

a = SI - 2(3+, 

( = SI + 0+ + V30-, 

7 = n + /3+ - V3/3-. (8) 

The momenta po, p± are canonically conjugate to fl, f3± respectively. For 
vacuum Bianchi IX or magnetic Bianchi VIo, the MSS potential V has the form 

V = c 2 e 2ba + e< + e 4 ~< - 2 (fle 2(a+( > + fl e 2(a+1 > + de 2 ^) . (9) 

Here a = X, b = 2, c = 1, and d = 1 for vacuum Bianchi Type IX, while 
a = 0, b = 1, c = v^Cj and d = — 1 for magnetic Bianchi Type VIo, and £ is a 
constant that depends on the strength of the magnetic field. In these models, 
the singularity occurs at r = — £1 = 00. From (||) and @, we see that V — > 
as t — ► 00 unless one of a, £, or 7 is w 0. When the potential itself is small 
(V w 0), (0) describes a "free particle" in MSS and is in fact approximately the 
Kasner solution. 

The standard algorithms for solving ODE's J2j| often employ an adaptive 
time step. The idea is to take large time steps where nothing much happens 
(e.g. in the Kasner regime) while taking shorter time steps when the forces are 
large (e.g. at a bounce). Unfortunately, in Mixmaster dynamics, the duration 
of the Kasner epochs increases exponentially in r as r — > 00 |^6| while the 
duration of the bounce itself is in some sense fixed p7| . This means that, 
although huge time steps may be taken in the Kasner segments, the time step 
must become very small at the bounces. Ideally, one would like the time step 
also to grow exponentially but this cannot be done with standard approaches. 
Thus, with a great deal of effort, standard methods can yield about 30 bounces 
for < t < 10 s in about an hour on a supercomputer. The application of the 
symplectic method to Mixmaster dynamics is straightforward. From (^), we 
have 

H K = -pl+p 2 + +p 2 _ ; H v = V(n,P + ,(3^. (10) 

With an adaptive step size and a 6th order version of (Eh, one can do slightly 
better (approximately a factor of three fewer steps) than 4th order Runge-Kutta. 
However, it is well known that a bounce off an exponential wall — the Bianchi II 
(or Taub 0) cosmology — is exactly solvable. If we first identify the dominant 
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Figure 1: Comparison of algorithms for Mixmastcr dynamics. A typical Mix- 
master bounce is shown in the anisotropy plane. Crosses indicates every 10th 
point on a 4th order Runge-Kutta evaluation of the trajectory, while circles 
indicate every 10th point for a 6th order standard symplectic evaluation. The 
filled squares indicate every point using the new algorithm. The inset shows the 
details closest to the bounce. Note that the new algorithm does not require a 
point at the apex of the bounce. 



exponential wall (say e ), then we find that the symplectic algorithm works 
for a different split of the Hamiltonian into two subhamiltonians [^0) . Let 

H = H 1 +H 2 (11) 

where (e.g. for Bianchi IX) 

H^-pl+pl+pi+e^-^ ■ H 2 =H v -e 4n - 8 ^. (12) 

Variation of H2 is exactly solvable as before, while variation of Hi yields equa- 
tions with solution 

Py (t + At) = Py {t) ; p_(t + At) =p-(t) ; 
y(t + At) = y{t) - 6 Py {t)At ; /3_(t + At) = /3_(t) + 2p_(t)At ; 
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Figure 2: A typical trajectory consisting of 268 Mixmaster bounces is shown 
projected onto the anisotropy plane. In terms of the rescaled variables /3+/|f2| 
and /?_/|f2|, the walls of the MSS potential and thus the bounce sites are at 
fixed locations. 



a(t + At) 



In 



— cosh4\/3(t + At -t ) 
E 



6p Q 



da 
~dt 



(13) 



where y = — 2ft + (3+ and E 2 = 3p 2 — p 2 _. This new symplectic algorithm 
provides an enormous advantage because the bounce is built in. The time step 
can grow exponentially. In a few minutes on a desktop computer, one can obtain, 
e.g., 268 bounces for 8 < t < 10 615 . Fig. [I] shows a comparison of the new and 
standard algorithms while Fig. |^ shows a typical trajectory. Since H\ is actually 
the exact Hamiltonian almost all the time, the accuracy of the method is much 
higher than the formal 6th order — we achieve machine precision. 

The BKL parameter u characterizes the angle in the anisotropy plane of the 
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Kasner epoch trajectory. Between the nth and n + 1st Kasner epochs, we have 

[ u n — 1 ; 2 < u n < oo 

U n +1 = < • (14) 

[ (u n - l)- 1 ; 1 < u n < 2 

All numerical studies || ^8|, p9[ p7j and a variety of analytic arguments |3(], Q 
have shown that one expects (|14|) to become ever more valid as r — > oo. With 
our algorithm in a double or quadruple precision code, one can evolve until 
deviations from ([TI]) disappear and then use the predicted values of u as a test 
of the accuracy of the code. For example [^0| , one discovers that the Hamiltonian 
constraint ((J?]) with H = 0) must be enforced although not necessarily at every 
time step. 

For our studies of spatially inhomogeneous cosmologies, it is important to 
keep in mind that 

(1) between bounces Mixmaster dynamics looks like the AVTD Kasner so- 
lution; 

(2) in a fixed time variable such as r, the time between bounces increases 
exponentially as r — * oo; 

(3) extraordinary accuracy can be achieved by symplectic methods when one 
of the subhamiltonians is dominant; 

(4) enormous gains in accuracy and computational speed can be made by 
using a "custom-designed" treatment of the bounce. 



4 The Gowdy Test Case 

As the simplest example of a spatially inhomogeneous cosmology, we consider 
the plane symmetric vacuum Gowdy universe on T 3 x R pl[ |32] | described by 
the metric 

ds 2 = e- x ' 2 e T / 2 {~e- 2T dT 2 + d9 2 ) 



[e p da 2 + 2e p QdadS + (e p Q 2 



')dS 2 



(15) 



where the background A and amplitudes P and Q of the + and x polarizations of 
gravitational waves are functions of r and < 8 < 2tt and periodic in 9. There 
is a curvature singularity at r = oo |32[ |[ ^3). The polarized case (Q = 0) 
has been shown to have an AVTD singularity B. The generic casejias been 
conjectured to be AVTD (except perhaps at a set of measure zero) 
we have verified to the extent possible numerically |2lJ, |35|, [37], |38|, |39, 
claim has been made that this model does not have an AVTD singularity 
which we believe to be incorrect E3. 

This model is especially attractive as a test case because Einstein's equations 
in our variables split into dynamical equations for the wave amplitudes (where 
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, a = d/da): 

P 1TT - e- 2T P, 90 - e 2P (Q 2 ~ e~ 21 -Q, 2 d ) = 0, (16) 

Q, TT -e- 2T Q,ee + 2 (P, T Q, T - e - 2T P, e Q, B ) = (17) 

while the Hamiltonian and momentum constraints become respectively first or- 
der equations for A: 

A, T - [P,? + e" 2T Pj + e 2P (Q 2 + e~ 2r Q 2 )] = 0, (18) 



X,e-2(P,gP, T +e 2P Q, 6 Q, T ) = 0. (19) 

Thus two problematical aspects of numerical relativity — preservation of the con- 
straints and solution of the initial value problem — become trivial |nj . The wave 
equations ( |l6] ) and ( fl7j ) may be obtained by variation of the Hamiltonian 



H = \ j ' d0 [ir 2 P 


\ J d9 [e~ 2r (P 2 + e 2P Q 2 )} = H K + H v (20) 



where irp and ttq are canonically conjugate to P and Q respectively. Variation 
of H K yields the AVTD solution 



P = In + v(t - T ) + ln[l + e - 2v{T - Ta) ] -> ut 



as t — > oo , 



1 e - 2 ' u ( T - T o) 

= Qo + ^ (l + e -2«(r-^)) "^o asr^oo, 
q _ e -2e(T-7n)\ 

TTp = V- —, — > V aST^OO, 

(1 + e- 2l, ( T - To )) 

TTQ = -2/*« (21) 

and is thus exactly solvable in terms of four functions of 6: fx, v > 0, Qo, and 
t . Hy is also (trivially) exactly solvable so that symplectic methods can be 
used @. 

If the singularity is AVTD, one would expect the spatial derivative terms to 
go to zero exponentially as r — > oo. However, if P — > vt, the term 

V 2 = e- 2 ^ 2P Q 2 (22) 
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Figure 3: P (solid curve), V\ (dash-dotted curve), and V2 (dashed curve) vs r at 
a fixed value of 9. Note that the slope of P, P, T , decreases after each interaction 
with V2 while P lT goes from negative to positive after each interaction with V\ . 
A continuation of this graph in r would show that V\ and V2 have permanently 
died off and that P continues to increase with fixed positive slope P, T < 1. 



in ( |l6| ) would grow rather than decay if v > 1. Thus Grubisic and Moncrief 
(GM) Q conjectured that, in a generic Gowdy model as t - +00, < v < 1 
except perhaps at isolated spatial points. If we consider generic initial data — 
e.g. P — 0, 7rp = v cos 9, Q = cos 9, ttq = — then we must ask how a generic 
Gowdy solution evolves toward the AVTD solution at each spatial point and how 
an initial P, T > 1 or P, r < is brought into the conjectured range. Typically, 
either V2 or 

V 1 =n 2 Q e- 2P (23) 
(where ttq = e 2P Q, T ) will dominate ( |l6| ) to yield either 

P, 2 T +n 2 Q e- 2P ~ K l (24) 

or 

Z 2 +Q,j e 2Z ^ k 2 2 (25) 

where Z = P — r. If P, T < 0, an interaction with V\ will occur to drive P, T 
to —P, T to yield P, T > 0. If P, T > 1, an interaction with V2 will occur to drive 
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(P, T — 1) to — (P, T — !)■ If this yields P, T < 0, a second interaction with V\ 
will occur, etc. When \P, T \ < 1, V2 disappears so that, after a possible final 
interaction with Vi, < P, T < 1 forever. A typical sequence of bounces is shown 
in Fig. 0. Fig. shows the maximum value of v on the spatial grid vs r. First 




1 10 100 



T 

Figure 4: Plot of f maa; vs t. The maximum value of v is found for two simula- 
tions with 3200 (solid line) and 20000 spatial grid points (circles) respectively. 
The horizontal line indicates v = 1. Continuation in r of the finer resolution 
simulation would show v max > 1 for a much longer r with v max < 1 eventually. 

we see that high values of r can be reached at which < v < 1 everywhere. 

Non-generic behavior can occur at isolated spatial points where either ttq or 
Q,g vanishes. In the former case, the absence of V\ where ttq — and its flatness 
where ttq w allow P, T and thus P to remain negative for a long time. Since 
Q, T = ttq e~ 2P , Q will grow exponentially in opposite directions on cither side 
of the points where ttq — producing a characteristic apparent discontinuity. 
On the other hand, if Q,g w 0, P, T can remain large for a long time causing 
a spiky feature in P. Both types of features sharpen and narrow with time. 
The features and their association with non-generic points are shown in Fig. ||. 
The presence of this non-generic behavior at isolated spatial points leads to a 
dependence of simulation results on the spatial resolution. The finer the spatial 
resolution, the closer will a grid point be to the non-generic point. Near these 
non-generic points, the generic process of approach to < P, T < 1 will occur 
but slowly since either ttq or Q,g w 0. The closer one is to a non-generic site, 
the longer this process will take. Thus a finer resolution code will have narrower 
spiky features at which it takes longer for P, T to move into the range [0, 1). Some 
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Figure 5: P (solid line) and Q (dashed line) vs 9 at t = 12.4 for the initial 
data set given here with v = 5 for < 6 < ir for a simulation containing 20000 
spatial grid points in the interval [0, 2ir]. The numbers on the graph refer to the 
most interesting features. Peaks #1, #2, and #3 in P are essentially the same 
in that they occur where Q,e ~ 0. Peak #4 shows an apparent discontinuity in 
Q where ttq w 0. 



evidence for this is seen in Fig. [| where the finer spatial resolution simulation 
diverges from the coarser one and will be considered in detail elsewhere p2| . 

Finally, we note that V% is analogous to the MSS potential. One may consider 
then application of the new Mixmaster algorithm to the Gowdy Hamiltonian 
20|) with H = H 1 +H 2 for 



Hi = I i> d9U 



2 , -2t+2P rt 2 



and 



'2 - 2 

These yield the exact solutions 

Q(r + Ar) = Q(r 



Q, 



H 2 = h<t> d8(7r 2 Q e~ 2P + e-^P*). 



(26) 
(27) 



ttq(t + At) 



Q(r),0 



tanh k(t — tq) 
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P(t + At) 



In 



Q(t), 6 



cosh k(t — tq) 



7Tp(r + Ar) = 1 — |«| tanhK(r — tq) 
from the variation of H\ (where £, k, and tq are functions of 9) and 



(28) 



P(t + At) 



P(r), 



ttq(t + At) = ttq(t), 



h p {t + At) = irp(T)+Tr 2 Q {T)e- 2P ^AT+±e- 2r (l - e - 2Ar ) P(r), 



mi 



Q(t + At) = Q(T)+ir Q (T)e- 2P ^A T 



(29) 



from the variation of i?2- Application of this algorithm is in progress. 
From the Gowdy test case, we learn that: 

(1) Since the singularity is AVTD, Hk dominates Hy asymptotically so our 
(current) algorithm is very accurate. 

(2) Non-linear terms in the wave equations act as potentials. In the Gowdy 
case, they drive the system to the AVTD regime as r — > oo with < V < 1, 
where the potentials permanently die out. 

(3) Non-generic points where ttq — or Q,g = lead to the growth of spiky 
features in P and Q. 

(4) Spiky features appear narrower with finer spatial resolution. 



5 U(l) Symmetric Cosmologies 

Given our understanding of the Gowdy model, we can move to spatially inho- 
mogeneous cosmologies with one Killing field rather than two, retaining a U(l) 
symmetry on T 3 x R ^]. These models can be described by five degrees of 
freedom {</?, to, A, z, x } and their respective conjugate momenta {p, r, p\, 
Pz, Px} which are functions of spatial variables u and v and time r. Einstein's 
equations may be obtained by variation of BlL E3| 



// = <j> <j) dudvTL 

' cfdudv(lp 2 z + \e^p\ + \p 2 + \e^r 2 - \p\ + 2p A ) 

- 2T I (f dudv {(e A e ab ) , ab - (eV b ) )0 A )6 +e A [(e~ 2z ) , u x, v - (e 
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+ 2eVV,a <p,t + \e A e'^e ab uo :a u, b } 
Hk + Hy = <j> <j> dudv Hk + j> j> dudvV. 



Here <p and u> are analogous to P and Q while A, x, z describe the metric e a b 
in the u-v plane perpendicular to the symmetry direction with 




This model is sufficiently generic that local Mixmaster dynamics is allowed. The 
U(l) Hamiltonian ( |30| ) has the standard symplectic form. Note that Hk consists 
of two Gowdy-like Hk 's and a free particle term so that Hk is exactly solvable. 
Hy is also (again trivially) exactly solvable although spatial differencing must 
be performed with care. In two spatial dimensions, there are a variety of ways 
to represent derivatives to a given order of accuracy. We currently use a scheme 
provided by Norton (43) to minimize the growth of short wavelength modes. 

Unlike the Gowdy model, the constraints and initial value problems must be 
considered. We find 

H° = H - 2 PA = (32) 

and 

H U = p z Z, u +PxX,u +PaA, u ~Pa,u +Pf,u +fW,a 

+±{[e 4z -(l+x) 2 ]p x -(l+x)p z } >v 

-±{[e 4 * + (l-x 2 )]p x -xp z }, u =0, (33) 

TL V = p z z, v +p x x, v +p A A, v -p\,v +P<P,v +rw, v 
-\{[e iz -{l-xf]p x + {l-x)p z }, a 

+| { [e 4z + (1 - x 2 )] Px - x Pz } , v = 0. (34) 

While a general solution to the initial value problem is not known, we use the 
particular solution obtained as follows: To solve the momentum constraints 
([H and set p x = p z = ip, a = uj, a =0 to leave p A A, a - p A , a = which may 
be satisfied by requiring p A — ce A . For sufficiently large c, the Hamiltonian 
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constraint may be solved algebraically for either p or r. In general, this leaves 
as free data the four functions x, z, A, and either r or p. Since there are four 
free functions at each spatial point, we expect generic behavior. 




Figure 6: Evolution of log 10 \V\ vs r for a polarized 17(1) cosmology at two 
representative spatial points. The solid lines are exponential fits. 



As a first case, we consider polarized J7(l) models obtained by setting 

r = u = 0. (35) 

This condition is preserved numerically as well as analytically. It has been con- 
jectured [H that polarized U(l) models are AVTD. This is reasonable because 
the Mixmaster potential-like term 

V Vu) =e~ 2T e A e ab e~^u;, a u J , b (36) 

is absent. Other spatial derivative terms decay as e A - 2r - 22 for the expected 
AVTD limits of the variables as r — > oo p3[ : 

Z — > ~V Z T , x — * x , p z — > -u z , 

Pa -»■ p° ) </? -?V T > -P > 

A -> A + (2-p°)r , p A ^p A , (37) 
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Figure 7: The maximum value of the change with time over the spatial grid 
of AVTD regime constants for a polarized U(l) symmetric model. Here v z = 
y / p 2 z /8 + e 4z p 2 x /2, VuJ = /8 + e 4 *Y 2 /2, c z = p z /2 + p x x, and c u = p/2 + rcu. 
Exponential decay is observed until the maximum values of the changes reach 
the level of machine precision. 



so that an AVTD singularity is consistent. 

Fig. U shows log 10 \ V\ vs r for typical spatial points. Thus we see the ex- 
pected exponetial decay. In Fig. ^, we see that the maximum values of the 
changes with time of the quantitites expected to be constant in an AVTD regime 
also decay exponentially as expected. These results will be discussed in detail 
elsewhere p3| . 

We emphasize here that the polarized U(l) models present no numerical 
difficulties. The absence of Vy w means that spiky features do not develop. 
The situation unfortunately changes for generic (unpolarized) 1) models. It 
appears that consistency arguments which suggest an AVTD singularity in po- 
larized U(l) models and restrict v to [0, 1] in the Gowdy models fail in generic 
U(l) models. This suggests that Vy w will always grow exponentially if the sys- 
tem tries to be AVTD producing a Mixmaster-like bounce. A bounce in the 
opposite direction will come, as in Gowdy models, from terms in H^. This 
bouncing could continue indefinitely. 

Unfortunately, the bounces off Vy w probably cause numerical instabilities 
that limit the duration of the simulations. Spatial averaging has been used 
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Figure 8: (a) log 10 V vs r for three representative spatial points of a generic 
U(l) symmetric cosmology, (b) The analogous quantity for a typical Mixmastcr 
evolution, (c) The analogous quantity for a Gowdy simulation. 



to improve stability but it is known to produce numerical artifacts. Typical 
evolutions of VVu> & t single spatial grid points are shown in Fig. || and compared 
to similar quantities in Mixmaster and Gowdy models and, of course, to Fig. ||. 
While generic U(l) models are clearly different from polarized ones, it is not 
clear yet whether the observed bounces are Mixmaster-like or will eventually die 
out as in Gowdy models. Recall that Mixmaster universes are very close to the 
AVTD Kasner solution between bounces. One is also concerned about numerical 
artifacts although it is not clear that standard tests are helpful. Fig. || shows 
two frames from a movie of V(u, v) vs r for two different spatial resolutions. 
Features are narrower on the finer grid. While this usually indicates artifacts, we 
recall that this is precisely what happens in Gowdy models where the resolution 
dependence is well understood. 



6 Future Directions 

In the search for the nature of the generic cosmological singularity, we have 
obtained convincing evidence that both the Gowdy universes and the polarized 
U(l) symmetric cosmologies have AVTD singularities. These results are found 
in simulations which present no numerical difficulties. In contrast, we cannot 
draw definite conclusions for generic U(l) models except to say that we have 
not found evidence that the singularities are AVTD everywhere. 

Presumably, numerical difficulties in generic U(l) models are due to Gowdy- 
like spiky features (which are less easy to represent accurately in two spatial 



dimensions). It is possible that the new Mixmaster algorithm 20 which can 
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Figure 9: Movie frames of V(u,v) for a generic U(l) symmetric cosmology for 
two fixed values of r. The upper frames precede the lower ones. The spatial 
coordinates u and v run from to 2tt in both directions. The frames on the 
left are from a simulation with 128 2 spatial grid points while the ones on the 
right have 256 2 spatial grid points. Both simulations have the same initial data. 
Note that there are features independent of spatial resolution as well as those 
that depend on the resolution. The background shade of gray indicates regions 
where V ~ 0. 



be adapted to the Gowdy model will help in the generic U(l) case where it can 
also be implemented. The hope is that better treatment of the bounces would 
give better local treatment of spiky features that arise from them. 

While we solve the constraints initially in U(l) models, we do nothing to 
preserve them thereafter. In the polarized case, they remain acceptably small 
(and in fact converge to zero with increasing spatial resolution). We have learned 
from the Mixmaster case that one must preserve the constraints [ p0[ . In fact, it 
is the kinetic part of the Hamiltonian constraint which restricts the exponential 
factor in VVw An error in the constraints could give the argument of the 
exponential the wrong sign leading to the observation of qualitatively wrong 
behavior. By starting closer to the singularity, we can supress some of the 
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numerical instabilities and study this controlling exponential. Studies of this 
type have provided some evidence that it is essential to solve the constraints. 
Work on implementing a constraint solver is in progress. 
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